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Abstract 

Rift Valley fever virus (RVFV) has been expanding its geographical distribution with important implications for both human and 
animal health. The emergence of Rift Valley fever (RVF) in the Middle East, and its continuing presence in many areas of Africa, 
has negatively impacted both medical and veterinary infrastructures and human morbidity, mortality, and economic endpoints. 
C\l Furthermore, worldwide attention should be directed towards the broader infection dynamics of RVFV, because suitable host, 
vector and environmental conditions for additional epidemics likely exist on other continents; including Asia, Europe and the 
Americas. We propose a new compartmentahzed model of RVF and the related ordinary differential equations to assess disease 
C/^ spread in both time and space; with the latter driven as a function of contact networks. Humans and livestock hosts and two 
species of vector mosquitoes are included in the model. The model is based on weighted contact networks, where nodes of the 
T— I networks represent geographical regions and the weights represent the level of contact between regional pairings for each set 

of species. The inclusion of human, animal, and vector movements among regions is new to RVF modeling. The movement 

of the infected individuals is not only treated as a possibility, but also an actuality that can be incorporated into the model. We 
have tested, calibrated, and evaluated the model using data from the recent 2010 RVF outbreak in South Africa as a case study; 
mapping the epidemic spread within and among three South African provinces. An extensive set of simulation results shows the 
O potential of the proposed approach for accurately modeling the RVF spreading process in additional regions of the world. The 
benefits of the proposed model are twofold: not only can the model differentiate the maximum number of infected individuals 
among different provinces, but also it can reproduce the different starting times of the outbreak in multiple locations. Finally, the 
exact value of the reproduction number is numerically computed and upper and lower bounds for the reproduction number are 
analytically derived in the case of homogeneous populations. 
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1. Introduction 



^— ^ Rift Valley fever (RVF) is a viral zoonosis with enormous health and economic impacts on domestic animals and humans 
. . [26], in countries where the disease is endemic and in others where sporadic epidemics and epizootics have occurred. An 
outbreak in South Africa in 1951 was estimated to have infected 20,000 people and killed 100,000 sheep and cattle [12, 35]. 
k> In Egypt in 1977, there were 18,000 human cases with 698 deaths resulting from the disease [12, 35]. While RVF is endemic 
5_i in Africa, it also represents a threat to Europe and Western hemispheres [7, 18]. In 1997 - 1998 Kenya experienced the largest 
^ recorded outbreak with 89, 000 human cases and 478 death. The first recorded outbreak outside of Africa occurred in the Arabian 
peninsula in 2000 - 2001 and caused 683 human cases and 95 deaths [17]. Tanzania and Somalia reported 1000 human cases 
and 300 deaths from an outbreak that was associated with above-normal rainfall in the region in 2006 - 2007 [17]. Rift Valley 
fever virus (RVFV) is generally distributed through regions of Eastern and Southern Africa where sheep and cattle are present 
[42]. It can cause morbidity (ranging from nondescript fever to meningo-encephalitis and hemorrhagic disease) and mortality 
(with case fatality rates of 0.2 - 5%) in humans [26]. The main economic losses of RVF in livestock arise due to abortion and 
mortality, which tends to be higher in young animals [9, 42], and bans on livestock exports during an epidemic [9, 4]. 

Rift Valley fever virus was first isolated from the blood of a newborn lamb in 1931 and later from the blood of adult sheep 
and cattle [44, 4]. Domestic ruminants and humans are among the mammalian hosts demonstrated to amplify RVFV [21] and 
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Aedes and Culex are believed to be the main arthropod vectors [7] . Rift Valley fever virus can be transferred vertically from 
females to their eggs in some species of the Aedes mosquitoes [18, 27]. The disease has been shown to be endemic in semi-arid 
zones, such as northern Senegal [45, 7, 28], and RVF epidemics often appears at 5 - 15 year cycles [28]. As noted earlier, RVFV 
has already spread outside Africa, to Yemen and Saudi Arabia [7, 6]. The species of vectors that are capable of transmitting 
RVFV have a wide global distribution [20] and there is therefore a distinct possibility for the virus to spread out of its currently 
expanding geographic range [9]. A pathways analysis [21] has shown that the RVF virus might be introduced into the United 
States in several different ways [21, 24] and that analysis identified several regions of the United States that are most susceptible 
to RVFV introduction. It is therefore desirable to develop effective models to better understand the potential dynamics of RVF in 
heretofore unaffected regions and then develop efficient mitigation strategies in case this virus appears in the Western hemisphere 
[18]. Such preparedness can help avoid a rapid spread of the virus throughout North America, as happened with the West Nile 
virus during the last decade [7, 18]. 

A RVF disease risk mapping model was developed by [1]. The authors observed sea surface temperature (SST) patterns, 
cloud cover, rainfall, and ecological indicators (primarily vegetation) via satellite data to evaluate different aspects of climate 
variability and their relationships to disease outbreaks in Africa and the Middle East [3,2]. The researchers successfully predicted 
areas where outbreaks of RVF in humans and animals were expected using climate data for the Horn of Africa from December 
2006 to May 2007. An ordinary differential equation (ODE) mathematical model was developed by [18]. The model is both 
an individual-based and deterministic model. The authors analyzed the stability of the model and tested the importance of the 
model parameters. However, neither human population parameters nor spatial (or, network) aspects are explicitly incorporated 
in the model. Another theoretical mathematical model on RVFV dynamic transmission was proposed [29]. This model is also an 
individual based model. The most important parameters to initial disease transmission and the endemic equilibrium have been 
carried out. 

In this paper, we present a novel model incorporating Aedes and Culex mosquito vector, and livestock and human host 
populations. Our model is based on weighted contact networks, where nodes of the networks represent geographical regions 
and weights represent the level of contact between regional pairs for each vector or host species. Environmental factors such as 
rainfall, temperature, wind and evaporation are incorporated into the model. For each subpopulation, a set of ordinary differential 
equations describes the dynamics of the population in a specific geographical location, and the transitions among the different 
compartments, after contracting the virus. We compute the lower and upper bounds of the reproduction number for homogeneous 
populations, explain their biological meaning, and numerically compare the bounds with exact values. 

We test, calibrate, and evaluate the model using the recent 2010 RVF outbreak in South Africa as a case study, mapping the 
epidemic spread in three South African provinces: Free State, Northern Cape, and Eastern Cape. An extensive set of simulation 
results shows the potential of the proposed approach to accurately describe the spatialtemporal evolution of RVF epidemics. 

The paper is organized as follows: 1) in Section 2, we describe our compartmentalized mathematical model, present the 
lower bound and upper bound of the reproduction number for homogeneous populations. 2) in Section 3, we introduce the case 
study using outbreak data from South Africa, 2010, 3) in Section 4, we conclude our work. In the Appendix, we show how we 
derive the bounds for the reproduction number for homogeneous populations. 

2. Compartmentalized Mathematical Model 

We have constructed Compartmentalized Mathematical Models based on the principle of RVFV transmission. The parame- 
ters used in the model are shown in Table 1 . 

2.1. Homogeneous Populations Model 

The principle of RVFV transmission between different species is shown in Figure 1 . All the Aedes ssp. and Culex spp. we are 
going to discuss only include the mosquitoes that are competent vectors of Rift Valley fever. In this paper, the Culex parameters 
are based on Culex Tarsalis mosquitoes and Aedes parameters are based on Ae. dorsalis mosquitoes. The main vectors, Aedes 
and Culex mosquitoes and the main hosts, livestock and humans are considered in the model. We use an SEI compartmental 
model in which individuals are either in a susceptible (S) state, an exposed (E) state, or an infected state (I) for both Aedes and 
Culex mosquitoes, and an SEIR compartmental model in which individuals are either in a susceptible (S) state, an exposed (E) 
state, an infected state (I), or a recovered (R) state for both livestock and human populations. Infectious Aedes mosquitoes can 
not only transmit RVFV to susceptible livestock and humans but also to their own eggs [18, 27]. Culex mosquitoes acquire the 
virus during blood meals on an infected animal and then amplify the transmission of RVFV through blood meals on livestock 
and humans [44]. Direct ruminant-to-human contact is the major (though not only) way for humans to acquire the infection 
[1, 11]. Accidental RVFV infections have been recorded in laboratory staff handling blood and tissue from infected animals 
[1]. Usually, humans are thought of as dead end hosts that do not contribute significantly to propagation of the epidemic [7]. 
There has been no direct human-to-human transmission of RVFV in field conditions recorded thus far [21]. The mosquitoes 
will not spontaneously recover once they become infectious [18]. Livestock and humans either perish from the infection or 
recover [18]. All four species have a specified incubation period [44]. The model is based on a daily time step. Aedes and Culex 
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mosquitoes are distributed among susceptible Sa, exposed Eg and infected 4 compartments. The subscript a — I denotes Aedes 
and a - 3 denotes Culex . The size of each aduh mosquito population is A^i = 5 1 + £1 + /i for adult Aedes mosquitoes and 
A^3 - S i + + It, for adult Culex mosquitoes. The livestock and human hosts contain susceptible Si,, exposed Ei,, infected I/, 
and recovered Ri, individuals. The subscript b - 2 and b - 4 denote livestock and humans respectively. The size of host popu- 
lations is Nh = S h+Eb+Ih+Rh- The four populations are modeled with a specified carrying capacity K\, K2, K3, K4 respectively. 
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Figure 1 : Flow diagram of RVFV transmission with each species, namely, Aedes mosquitoes, Culex mosquitoes, livestock, and humans homogeneously mixed 
(the solid lines represent transition between compartments and the dash lines represent the transmission between dilferent species) 



2.1.1. Aedes Mosquito Population Model 
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Parameter 


Description 


Value 


Units 


Source 


Pn 


contact rate: Aedes to livestock 


0.002 


1/day 


[18] 


Pii 


contact rate: livestock to Aedes 


0.0021 


1/day 


[18] 




contact rate: livestock to Culex 


0.000003 


1/day 


[18] 


P-ii 


contact rate: Culex to livestock 


0.00001 


1/day 


[18] 


P\A 


contact rate: Aedes to humans 


0.000046 


1/day 


Assume 


PlA 


contact rate: livestock to humans 


0.00017 


1/day 


[29] 




contact rate: Culex to humans 


0.0000001 


1/day 


Assume 


J2 


recover rate in livestock 


0.14 


1/day 


[18] 


74 


recover rate in humans 


0.14 


1/day 


[39, 38, 37] 


d\ 


death rate of Aedes mosquitoes 


0.025 


1/day 


[18] 


d2 


death rate of livestock 


1/3650 


1/day 


[18] 




death rate of Culex mosquitoes 


0.025 


1/day 


[18] 


d4 


death rate of humans 


1/18615 


1/day 


[39, 38, 37] 


bi 


number of Aedes eggs laid per day 


0.05 


1/day 


[18] 


bi 


daily birthrate of livestock 


0.0028 


1/day 


[18] 


b3 


number of Culex eggs laid per day 


weather dependent 


1/day 


[19] 


bA 


daily birthrate of humans 


1/14600 


1/day 


[39, 38, 37] 




incubation period in Aedes mosquitoes 


6 


days 


[18] 


1/62 


incubation period in livestock 


4 


days 


[18] 


IM 


incubation period in Culex mosquitoes 


6 


days 


[18] 


1/64 


incubation period in humans 


4 


days 


[44] 




mortality rate in livestock 


0.0312 


1/day 


[18] 




mortality rate in humans 


0.0001 


1/day 


[39, 38, 37] 




transovarial transmission rate in Aedes 


0.05 




[18] 


1/^?1 


development time of Aedes 


15 


days 


[18] 


03 


development rate of Culex 


weather dependent 


1/day 


[19] 




carrying capacity of Aedes mosquitoes 


1000000000 




[32] 


Ki 


carrying capacity of livestock 


10000000 




Assume 


K, 


carrying capacity of Culex mosquitoes 


1000000000 




[32] 


Ka 


carrying capacity of humans 


10000000 




Assume 


f 


fraction of those working with animals 


0.82 




[31] 


T 


return rate 


3 


times/day 


[5] 


P 


reduction in cj^j due to infection 


1 

2 




Assume 



Table 1 : Parameters of the compartmentalized mathematical model 
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^ ^0,{Pi + Q,)-d,NiNdKi (6) 
at 

where: 

f i=the number of uninfected Aet/ei mosquito eggs 
!2i=the number of infected Aedes mosquito eggs 
S 1 =the number of susceptible Aedes mosquitoes 
£'i=the number of exposed Aedes mosquitoes 
Ii =the number of infected Aedes mosquitoes 
A^i =the total number of Aedes mosquitoes 

The above model is a modified SET model with compartments P and Q. Compartments P and Q represent uninfected eggs 
and infected eggs respectively. The total number of eggs laid each day is b\N\ with biq\I\ infected eggs and b\N\ - b\q\I\ unin- 
fected eggs [18]. After development period, 6\Pi develop into susceptible adult mosquitoes and 9i Q\ develop into infected adult 
mosquitoes [18]. There are d\X\N\IK\ mosquitoes removed from compartment X due to natural death. Compartment X can be 
P, Q, S, E, and I here. The number of Aedes mosquitoes infected by livestock is denoted by I32\S \I2IN2 which is proportional to 
the density of infected Aec/e* mosquitoes [18]. After incubation period, e\E\ Aedes mosquitoes transfer to infected compartment 
[18]. 

2.7.2. Culex Mosquito Population Model 



: biNi - 6iPi (7) 

■ O3P3 -/323S3I2/N2 - d3S3N3/K3 (8) 

■ /323S3I2/N2 - S3E3 - d3E3N3/K3 (9) 

■ E3E3-d3l3N3/K3 (10) 
■■ 63P3-d3N3N3/K3 (11) 



dP3 
"dT 
dS3 
"dT 

dE3 

~dt~ 

d/3 

dT 
dA^B 
"dT 

where: 

f 3=the number of uninfected Culex mosquito eggs 
5'3=the number of susceptible Culex mosquitoes 
£'3=the number of exposed Culex mosquitoes 
/3=the number of infected Culex mosquitoes 
Af3=the total number of Culex mosquitoes 

Besides compartment S , E, I, compartment P is added to represent uninfected eggs. Only uninfected eggs are included be- 
cause the female Culex mosquitoes do not transmit RVFV vertically [18]. The total number of eggs laid each day is b3N3. There 
are d3X3N3/K3 Culex mosquitoes removed due to natural death. Compartment X can be P, S, E, and I here. After development 
period, 63P3 eggs develop into susceptible adult Culex mosquitoes and become secondary vectors [18]. The number of infected 
Culex mosquitoes from contact with livestock is denoted by P23S3I2IN2 which is proportional to the percentage of infected live- 
stock [18]. After incubation period, E3E3 Culex mosquitoes transfer from exposed compartment to infected compartment [18]. 

2.1.3. Livestock Population Model 



At 
d£2 
At 
Ah 
At 
AR2 
At 



-■ b2N2 - d2S2N2/K2 -/312S2I1/N, -/332S2I3/N3 (12) 

-■l3nS2hlNi +/332S2I3/N3 - £2^2 - d2E2N2/K2 (13) 

: e2E2 - y2h - f^ih - 1/2^2-^2/-'^2 (14) 

-.y2h-d2R2N2lK2 (15) 
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— ^ = b2N2 - d^NiN^IKt - ii2h 



(16) 



where: 

52=the number of susceptible livestock 
£'2=the number of exposed livestock 
/2=the number of infected livestock 
A^2=the total number of livestock 



There are h2N2 livestock born, d2^2NilK.2 livestock removed due to natural death [18], and //2^2 livestock dying of the infec- 
tion each day [18]. Compartment X can be S, E, I, and R here. Following incubation period, £2^2 livestock transfer from exposed 
compartment to infected compartment [18]. The number of livestock infected by Aedes mosquitoes and Culex mosquitoes are 
denoted as PnSihlNi and f332S2h/Ni respectively [18]. Following infection period, 72^2 livestock recover from RVFV infec- 
tion [18]. 



2.7.4. Human Population Model 



dS4 
"dT 
d£4 
"dT 
d/4 
dT 
di?4 
"dT 
dA^4 
"dT 



■■ bANA-PuSAhlNi - fl32ASAhlN2 - ^348 Ah/Ni - d^SANAlK^ 
■■ PxaSaIxIN, + f/324S4h/N2 +/^34S4h/N3 - af4£4A^4/^4 - £4^4 

■ e4£'4 - 74/4 - yU4/4 - d4l4N4/K4 

■ 74/4 - d4R4N4/K4 

■ b4N4 - d4N4N4/K4 - JU4/4 



where: 

54=the number of susceptible humans 
£'4=the number of exposed humans 
/4=the number of infected humans 
N4-the total number of humans 



(17) 
(18) 
(19) 
(20) 
(21) 



There are b4N4 humans born, d4X4N4/K4 humans removed from compartment X due to natural death, and JU4/4 humans 
dying of RVFV infection each day. Compartment X can be S, E, I, and R here. The number of humans that acquire the infection 
from Aedes mosquitoes, Culex mosquitoes, and livestock is /314S4I1/N1, /334S4I3/N3, and f/324S 4I2/N2 respectively. We assume 
only those who work with animals can be infected by animals. Therefore, a coefficient / (0 < / < 1) which represents the 
fraction of humans working with animals is multiplied by j6245' 4/2/^2- After incubation period, e4£'4 humans transfer to infected 
compartment and 74/4 humans transfer to recovered compartment after infection period. 



2.1.5. Environmental Parameters for Culex 

The equation (22) is used to model the development rate of Culex mosquitoes [19]. The daily egg laying rate expressed in 
equation (23) is a function of moisture [19]. Moisture in equation (24) is obtained by summing the difference of precipitation 
[30] and evaporation (mm) [25] over the proceeding 7 days [19]. In the equations (22) to (25), A, HA, HH, K, TH, E,na.x, Eyar, 
Emean, t>o are parameters [19] which are described in Table 2. This model is specific for West Nile virus model in 2010 in the 
northern US. 



e3(Temp, t)^A. (Tempit) + K) ^ ^^/^[ i^^v*( 298^5 - r.JgH^)] 
298.15 i^exp[J%Hj^-r^K)^ 
Emax 

b3(Temp, precipitation, t) ^ bo + — Mois,urei,)-E„,ean -. ^23) 

' ^ Evar -I 
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Parameter 


Description 


Value 


Source 


A 


parameter of the development rate 


0.25 


[19] 


HA 


parameter of the development rate 


28094 


[19] 


HH 


parameter of the development rate 


35692 


[19] 


TH 


parameter of the development rate 


298.6 


[19] 


bo 


minimum constant fecundity rate 


3 


[19] 


^max 


maximum daily egg laying rate 


3 


[19] 




value at which moisture index=0.5£'„,a j 





[19] 


p 


the variance of the daily egg laying rate 


12 


[19] 



Table 2: Parameters of the model for Culex 



Moisture{t) — ^ participation{D) — evaporation(D) (24) 

D=t-6 

100(Temp(t) + 0.006/!)/(100 - latitude) 



Evaporation(t) 



80 - Tempit^ 
\5{Temp{t) - TAt)) 



80 - Temp(t) 

Where: 

Temp{t) =air temperature in units of "C [25] 

latitude =the latitude (degrees) [25] 

Td(t) =the mean dew-point in units of "C [25] 

h -the elevation (meters) [25] 

K - Kelvin parameter [25] 



mm /day (25) 



2.7.6. The Reproduction Number for Homogeneous Populations 

The reproduction number Rq is defined as: " The average number of secondary cases arising from an average primary case 
in an entirely susceptible population" [14]. The reproduction number is used to predict whether the epidemic will spread or die 
out. There are several methods used to compute Ro- One of these methods computes the reproduction number as the spectral 
radius [14? ] of the next generation matrix [14? ]. 

The next generation matrix is defined as and the matrices F and V are determined as: 

F = [—^ ], V = [— ] 

oxj oxj 

where xj is the number or proportion of infected individuals in compartment j, j - 1, 2, 3, ... , m, m being the total number of 
infected compartments, xq is the disease free equilibrium vector, is the rate of appearance of new infections in compartment 
and fi - 'Yr - ^+ with i^r denoting the transfer of individuals out of compartment / and "f.^ denoting the rate of transfer 
of individuals into compartment / [16]. The {i, j) entry of F represents the rate at which infected individuals in compartment j 
produce infected individuals in compartment / [16]. The {j, k) entry of V ' represents the average time that an individual spends 
in compartment j, where /, j,k - 1,2, 3, ... , m [16]. Finally, the (/, k) entry of FV^^ represents the expected number of infected 
individuals in compartment / produced by the infected individuals in compartment ^ [16]. 

For our homogeneous population model, we found that 



Ri; ^Ro^R^ + qi (26) 

where 



^0 - 



£2 



;[ 



^ (bi + S2)(b2 + 72 + mVbiibi +si) b^{bj. + e^) 



(27) 



See the Appendix for the derivation details, the biological interpretation, and the comparison among exact values and bounds 
for the reproduction number. 
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2.2. Meta-Population Model 

A meta-population model is a model with several subpopulations. It assumes homogeneity within each subpopulation and 
heterogeneity among different subpopulations. The Aedes and Culex mosquitoes in location i (i - 1, 2, 3, ■ ■ ■ , n), are distributed 
among susceptible Sai, exposed E^i and infected /«, compartments. The subscript a - I denotes A eiie.s and a = 3 denotes Culex . 
The size of each adult mosquito population in location / is A^i, = 5 ii+Zsii+Zi, for adult Aet/ei mosquitoes andA^3, = S^i + E^i + hi 
for adult Culex mosquitoes. The livestock and human hosts contain susceptible S hi, exposed Ehi, infected 7^, and recovered Rt,i 
individuals. The subscript b - 2 and b - 4 denote the livestock and humans, respectively. The size of host populations in 
location ; is N2i - 821 + £2/ + hi + Rn for livestock hosts and Nm - Sm + E41 + + Rm for human hosts. The four populations 
are modeled with a specified carrying capacity K\, K2, K^, K4 respectively. 

2.2.7. Movement between Nodes 

We used weighted networks for each compartment of the four species as is shown in Figure 2. The superscripts of oj on the 
left hand side of equations (28), (29), and (31) represent the movement of different species. The number '1' in the superscript 
means the movement of Aedes or Culex population, '2' means the livestock movement, and '3' means the human movement in 
the superscript. The subscript ij of aj,y means that the direction of the movement is from location / to location j. The difference 
in the thickness of the lines represent the difference in weight. Thicker lines represent the larger weight. The weight for each 
population is between and 1. RVFV has been documented to be spread by wind [35]. Wind dispersal of mosquitoes has 
changed geographic distribution and accelerated the spread of RVFV to new geographic areas [21]. Some locations can become 
secondary epidemic sites after the virus has been introduced (especially in irrigated areas, e.g. Gazeera in Sudan or rice valleys in 
the center of Madagascar) [28]. Livestock trade and transport also can affect the geographic distribution of RVF [7]. One critical 
objective in developing effective models is to determine the major factors involved in the disease spreading process. Therefore, 
we parameterize the weight due to mosquito movement with wind [21, 8], livestock movement due to transportation to feedlots 
or trade centers [40], and human mobility due to commuting [5] as shown in equations (28), (29), and (31), respectively. The 
movement rate of infected livestock is reduced due to infection [44]. We use the wind data [41] in Bloemfontein, which is the 
capital of Free State, as the wind of Free State Province, that of Kimberley, which is the capital of Northern Cape, as the wind of 
Northern Cape Province and that of Grahams town, which is the center of Eastern Cape Province, as the wind of Eastern Cape 
Province. The distance vector is calculated with longitude and latitude in the center of each location. The number of animals sold 
[37] and the number of livestock in the feedlots [34] are factors of weight for livestock movement. Distance, human population, 
commuting rate, and return rate [39] affect the weight for human movement. Weight for mosquito movement is decided by 
distance and the projection of wind in the direction of distance vector [8]. 

4 

Here: 

Wi= the wind vector in location / [8] 
Dij = the distance vector from location / to location j 
<x)\j{t) - the weight for mosquitoes moving from location / to location j 
cojj(t) - the weight for livestock moving from location / to location j 
(Tij(t) = the number of commuters between location / and location j 
FMi - the number of animals in markets and feedlots in location / 

2.2.2. Aedes Movement between Nodes 

^ =bi{Nu-qiIu)-OiPu (33) 



ci- 



C2 



- C3- 



WrDjj 1 

\Dij\ \Dij\ 
FMj 1 

FMi IT?;-! 



N4i 



- Z 



(28) 
(29) 

(30) 
(31) 
(32) 
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(a) Mosquito movement network (Mosquitoes can 
move from node (' to node ji, ji, and 73 and vice versa 
due to wind. We assume mosquitoes do not return to 
the node tliey are from.) 




(b) Livestock movement network (Livestock can move 
from node (' to node j\, j2, and 73 and vice versa due 
to trade. We assume livestock do not return to the node 
they are from.) 




(c) Human movement network (Humans can commute 
from node i to node 71, j2, and 73 and vice versa. We 
assume humans return to the node they are from.) 

Figure 2: Network graphs with node which has tliree neighbors as an example 
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dt 

dSu 
dt 

dEu 
dt 

diu 

dt 

dNu 
dt 



= biqJii - 0iQu 

n n 

^BiPu-l3iiSuhilN2i-diSuNulK,+ J] oj],S ij - J] oj]jS u 

ti n 

^ PixS uhilN2i- e,Eu- diEuNuIKi + ^ - Yj 

n n 

= 0iQu + siEii-diIuNu/Ki+ ^ 2 ojllu 

n n n 

j=i,;5'' ;=ijVi j=i,>5'' 



j=lj*i j=lj*' 

The change in the number of Aedes mosquitoes due to mobility in compartment X is given as j^i 
[22]. 

2.2.3. Culex Movement between Nodes 



dt 

dt 

dEjj 
dt 

div 
dt 

dNii 
dt 



n n 

-■OiPii- 132^8 iihilN2i-diSiiNiilKi+ ^ oj]iS^j- J] w^/s/ 

n n 

-■l323SiihilN2i-SiEii-diEiiNiilKi+ ^ oj]fiij - ^ ^.£3/ 

n n 

-■ eiEj,i-d^hiN3ilKi+ Y ^ i^\jhi 

n n n 

-.03P3i-d3N3iN3i/Ki+ Y <^)iS3j- Y ^,7^3/+ Y "^'ji^^j 



j=Uj*i j=l,j*' j=lj*i 



j=l,j*i j=^,i*' j=l-j*i 

The change in the number of Culex mosquitoes in compartment X due to movement is given as 2"=i ^)/^3j 
[22]. 

2.2.4. Livestock Movement between Nodes 



dS " " 

= b2N2i -pnS2ihi/Nu -P32S2ihi/N3i - d2S2iN2i/K2 + Z 4.52) - Z colS2i 



dE2i 
dt 

dhi 
dt 

dR2i 
dt 



-- Pi2S2ihilNu+ l332S2ihilN3i - e2E2i- d2E2iN2ilK2+ Yj '^lEij- X 4-^2' 

n n 

-■ p Yj ^ji^y - P Yj ^^z'^^' ~ '^2hiN2i/K2 + E2E2i - y2hi - P-ihi 

n n 
-. ^iRy- Yj '^l'^2i+nl2i-d2R2iN2i/K2 
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AM n n n 

^b2N2i-d2N2iN2ilK2-H2hi+ J] W^.^j; - ^ OJ^jS 2i + Yj "^l^^j 

j=^:i*' .i='^J*i j=lJ*' 

n n n n n 

j=lj*' j=l.j^i j=lji=i j='^,j*' j=l,j^i 

The change in the number of animals due to movement in susceptible, exposed, and recovered compartment is Y/j^i ji^i <^^jp^2j~ 
ji^i '^fj^ii [22] for livestock. Concerning the animals in the infected compartment, we assume that the movement rate of the 
infected livestock is p (0 < p < 1) of hvestock in other compartments. This value of the movement rate has been selected in the 
absence of further information. 

2.2.5. Human Movement between Nodes 



—— = b^NAi - d^SAiNM K4 ; 

d/ 1 + cr,/T 1 + cr,/T 1 + cr,/T 

1 + 0-,/t .4-^. 1 + cr,/T 4-^. 1 + cr,/T ^ ' 

d£4/ ^ PxASAihilNxi ^ j324fS4ihj/N2i ^ PuSAihilN^ ^ y PxaS AihjlNxjtTjjlT 

At 1 + (T,/t 1 + cr,/T 1 + cr,/r .4^^^. 1 + cr,/T 

Z" /32AfS4il2i/N2iO-ii/T B^aS Aihi/NiiO-ii/T 

, ' ' + y ; ' ' -dAEAiNAi/KA-SAEAi (50) 

.... 1 + CTij T . 1 + O-j T 

^ = BAEAi - jAhi - HAhi - dAhiNAi/KA (5 1) 

at 

^ = jAhi - dARAiNAi/KA (52) 
df 

— ^ = bANAi - dANAiNAilKA - fiAhi (53) 
at 

The humans from location / can stay in location / or move to location j at time f [5]. The number of humans infected by 
Aedes mosquitoes, Culex mosquitoes and livestock is PiAiSAuj^ + j^jSAijj^) [5], /33a(S au]^ + liy=i ji^i Ai j j^) [5], and 
fl32A{SAu^^ + Ij"j=i,j^iSAij]^^) [5] respectively. 

where: 

54,,= the number of humans that are from location / and stay in location / at time t [5]. 
Saij- the number of humans that are from location / and stay in location j at time f [5]. 
(o]j= the commuting rate between subpopulation / and each of its neighbor j [5] 
Wj - daily total rate of commuting for population / [5] 

The change in the number of susceptible humans that are from location / and stay in location / is given [5] by the following 
expression. 

dS AH 

dt 



Y ^'^4,7 " Y 4/4,7 



The change in the number of susceptible humans that are from location / and stay in location j is given [5] by the following 
expression. 



dS 



^'V _ , .3 



- W; :S Aii - tS Aij 



dt 

We can get the solution of 54,-, and Saij through the above two equations at the equilibrium. 



Sau = , , (54) 

1 + (jJi/T 

Sau = 1 , 4/'^ (55) 
1 + a>i/T ■' 
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3. Case Study: South Africa 2010 

We have used data from the South African RVF outbreak in 2010 as a case study. 

3.1. Incidence Data Analysis 

Outbreak data for animals are obtained from [15, 43], while outbreak data for human subpopulations are collected from 
[13, 31]. As far as animal data is concerned, we chose to analyze RVF incidence in the sheep population. Because the granularity 
of human incidence data is provided at Province level, each node in the network represents a province. We selected three 
provinces: Free State (location 1), Northern Cape (location 2) and Eastern Cape (location 3), because they had the highest levels 
of RVF incidence for humans. The curves of the incidence data are shown in Figure 3 using green histograms, while the red 
curves represent simulations obtained with our model. From the data in Figure 3, it is possible to observe that the epidemic 
started first in the Free State Province and later in Northern Cape Province. The sustained heavy rainfall likely triggered the 
outbreak, causing infected eggs to hatch in the Free State Province. Additionally, the number of animal and human cases in 
Eastern Cape Province is smaller than the other two provinces. 

3.2. Sensitivity Analysis 

The three parameters ci, C2 and C3 are estimated using the least square approach, which is based on minimization of errors 
between the incidence data of humans and the percentage of humans calculated by the mathematical model. At first, we establish 
an objective function. At each sample time, we calculate the difference between the number of humans calculated by differential 
equations and that reported [36] during outbreaks in three provinces of South Africa from January, 2010. We calculate the square 
of each difference. Then, we add all the squares for each location in each day together to obtain the objective function as is shown 
below. Minimization of the objective function is initiated by providing initial values cio, C20 and c^o for each parameter. The 
differential equations are solved with each set of the parameters and the square errors between the number of infected humans 
obtained from the objective function and those from incidence data are calculated. The parameters ci, C2 and C3 we used in the 
model are ci = 0.009, C2 = 0.05 and C3 = 0.005. 



In the equations above, 
n - the number of nodes 
to =starting time 
tf =end time 

hiit) =human prevalence calculated by the model 
PRiiiit) =human prevalence reported 

To conduct a sensitivity analysis of the parameters c\, C2, and C3 in equations (28), (29), and (31), we have changed each 
parameter within +10% of the values ci - 0.009, C2 - 0.05 and C3 = 0.005, keeping the other parameters constant. This analysis 
allows an evaluation of the impact of uncertainties in the parameter estimations. The percentage of infected humans obtained 
from simulation with ci = 0.009, C2 - 0.05 and C3 = 0.005 is represented as loiit). The percentage of infected human obtained 
from simulation with the parameters within ±10% bound is represented as /4,(f)- The relative errors between the fractions of 
infected humans are calculated for each set of parameters, in each location, at time f as \ The relative errors and the 

lower bound and upper bound of the parameters are shown in Figure 4. 

All the values of relative errors shown in Figure 4 are smaller than 10%, proving the model robustness with respect to limited 
uncertainties in the parameter estimation. The rest of the parameters such as contact rate /3i2, Pix^ Pn^ Pn, death rate di, dy, 
and recovery rate 72 are the most significant parameters in [18]. Similarly, /3i4, ^624, j034 and 74 are also the most significant 
parameters in this model. 

3.3. Analysis of Simulation Results 

To explore the behavior of RVFV, we conducted numerical simulations of an open system considering movement of the 
four species among different locations. To test the validity of the model, we changed some parameters in the weights to see 
the impact of each variation. If the number of infected eggs 2ii - 10, Q12 = 0, and Qn - 0, at the beginning infected eggs 
only exist in location 1, Free State. However, our model considers movement of mosquitoes to other locations with wind. As a 
consequence, infected animals and humans appear in all three locations as is shown in Figure 5. Therefore, the infection spreads 
due to movement of the four populations. 
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(56) 
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gure 3: Simulation results and incidence data from January, 2010 in South Africa (bars are data and lines are simulation results) 
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(a) Relative error of the number of (b) Relative error of the number of (c) Relative error of the number of 
infected humans with different value infected humans with different value infected humans with different value 
ofci ofc2 ofc3 



Figure 4: The relative eiTor of the number of infected humans with changing one of the parameters c\,C2, and C3 




Time(days) Time(days) Time(days) 



(a) Free State province (b) Northern Cape province (c) Eastern Cape province 

Figure 5: Simulation results with nonzero movement weights (the solid line represents livestock with y-axis on the left and the dash line represents humans with 
y-axis on the right) 



If we also assume that at the beginning infected eggs only exist in location 1, Qn = 10, Q\2 - and = 0, and movements 
of the four species from one location to another are not allowed, ci = 0, C2 = 0, and C3 = 0, then infected animals and humans 
will not appear in location 2 and location 3 as is shown in Figure 6. We can test the mitigation strategy of movement ban with 
this model. We performed the simulations to reproduce the RVF outbreak in the three South African Provinces. The simulation 
results and the incidence data are shown in Figure 3. The model can differentiate the maximum number of infected individuals 
among the three different provinces, it also reproduces the different starting time of the outbreak in the three locations. With a 
homogeneous population model, such as the one in [18], the spatial differentiation is not possible. 

The animal incidence curves provided by the model were always an overestimation of the data, since underreporting is 
very common during outbreaks. Finally, our approach in which the fractions of each subpopulation in each compartment are 
expressed as continuous variables, requires a large number of cases to be accurate. For this reason, the incidence data for location 
3, Eastern Cape, are better approximated by a stochastic model. The model has shown the ability of fitting the data. The starting 
time and trend of outbreak dynamics have been reproduced by the model. 



4. Conclusions 

A meta-population, network-based, deterministic RVF model is presented here. The animal, human and mosquito movement 
and their spatial distribution are considered by the model. The model successfully describes a real outbreak dynamics of RVFV, 
taking into account space and movement. When considering n locations or nodes (n >= 1), there are 21 differential equations 
and 21 *« variables in our model, while there are only 14 equations with 14 variables in the model presented in [29] and 16 equa- 
tions with 16 variables presented in [18]. Greater accuracy of our model is obtained at the cost of an increased complexity. The 
novelty of our model is that it considers a weighted contact network to represent the movement of four species. Subpopulations at 
the node level are also incorporated in our model. Additionally, parameters representing mosquito propagation and development 
are not constant but are the functions of climate factors. The model has been evaluated using data from the recent outbreak in 
South Africa. We reproduced not only the starting time but also the trend of RVFV transmission with time in different locations. 
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Figure 6: Simulation results with cl = c2 = c 3 = (the solid line represents livestock with y-axis on the left and the dash line represents humans with y-axis on 
the right) 



The model has shown to be very promising notwithstanding the limitation of the data. Due to the flexibility and accuracy of 
the proposed model, we can test and design multiple and different mitigation strategies in diff'erent locations at different times. 
The lower bound and upper bound of the reproduction number for homogeneous populations are shown to be very close to the 
exact value, and they provide insights on the biology of the spreading process. Future work in follow-up mathematical models 
includes the development of a stochastic model, the study of the impact of climate changes on the epidemiology and control of 
RVF, and the improvement of the mosquito movement model considering diffusion equations. Moreover, the carrying capacities 
of mosquitoes will be considered dependent on climate factors in the future. 
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Appendix 

Exact Computation ofRo 

We compute Rq as the spectral radius of the next generation matrix of the entire system [14, 29], Rq - p{FV^^). Before 
applying the method we need to verify that the five assumptions in [16] are satisfied [23]. First, the equations in the system are 



reordered so that the first m (m - 9) compartments correspond to infected individuals. 

- biqJi - 0iQi (57) 

^ ^PiiSih/Ni-siEi-diEiNi/Ki (58) 

^^6iQ,+siEi-dihNJKi (59) 
at 

^ = y8i25 2/1 /Ni + /532S 2/3 /A^3 - S2E2 - d2E2N2/K2 (60) 
at 

d/2 

— = S2E2 - 72^2 - J"2^2 - d2hN2lK2 (61) 

at 

^ = l323S3h/N2 - £3^3 " d.EiNj/K, (62) 
df 

^^s^E^-d.hNi/K, (63) 
at 

^ =l3i4S4h/m +f/524S4h/N2+l334S4h/N,-d4E4N4/K4-S4E4 (64) 

df 



15 



d/4 
dT 

dt 
dPi 

"dT 
dSi 
"dT 
dS2 
"dT 
d53 
"dT 
dS4 
"dT 
d/?2 
"dT 
d/?4 

"dT 



b^N^ - 03P3 

9iPi -/32iSil2/N2-diSiNJKi 

b2N2 - d2S2N2/K2-/3nS2li/Ni -Pi2S2hlN3 

6*3^3 -l32iSil2lN2 - diSiNilK^ 

b^NA -PuSAhlNi - ,fl32ASAhlN2 -/33ASAh/Ni - dASANA/KA 

y2h - d2RiNilK.2 

74/4 - dA^A^Al^A 



(65) 
(66) 
(67) 
(68) 
(69) 
(70) 
(71) 
(72) 
(73) 



The above system can be written as /^(x) - -^kix) - %(x), k- 17 
where 

x^[Qi El Ii E2 h £3 h Ea h Pi 
is the number of individuals in each compartment. 

Xs- [6? /? £° q £3" ^3 ^4 ^4 P' 

= [0 00000000 



Si S2 S3 Sa Ri RaY , 



di0i 



pO 

fc|£3 

1/363 



5? 



Si 

b2K2 

di 



S'3 

d. 



cO 
'^4 

d^ 



Rl RlY 
0]^ 



is the set of disease free states. 

^(x), y^ix), and y^(x) are given in the following. 



^1 = 




rf = 


01 gl 






^2 = 


I321S1I2/N2, 


■^2^ = 


fii^i +diEiNilKi 






^3 = 


0, 




dihNiIKi, 






=^4 = 


j3l2S2h/Ni +/332S2h/N3, 




S2E2 + d2E2N2/K2 






=^5 = 


0, 




72^2 + f^2h + d2l2N2/K2, 




= e2£'2 


=^6 = 


j02353/2/A^2, 




S3 E3 + d3 E3N3 1 K3 






=^7 = 


0, 




d3l3N3/K3 




= e3£'3, 


=^8 = 


PuSaIiINi +fl32ASAhlN2, 
+P34SAh/N3, 




dAEANA/KA + saEa 






^9 = 


0, 




jaU + JU4/4 + dAhNA/KA, 




= e4£'4 


=^10 = 


0, 


-^10 




'^10 


^biNi 


=^11 = 


0, 


■^n = 


03P3, 


'^ll 


= /73A^3 


=^12 = 


0, 


-^12 = 


f^2iSil2/N2 + diSiNi/Ki, 


-^12 


= 0lPi 


=^13 = 


0, 


-^13 


d2S2N2/K2 +/3nS2ll/Ni +P32S2h/N3, 


^^13 


= b2N2 


=^14 = 


0, 


-^14 


j323S3l2/N2 + d3S3N3/K3, 


-^14 


^03P3 


=^15 


= 0, 




^l3iASAhlNi+ fli2ASAhlN2, 
+/33ASAh/N3 + dASANA/KA 




= /74A^4 


=^16 = 


0, 


-^16 


d2R2N2/K2 




= 72^2 


=^17 = 


0, 




dARANA/KA, 


'^17 


= 74^4 



As it can been easily seen, the following five assumptions [16] are satisfied. 
(Al) if X ^ 0, then .^i, y.^ , ^Ofoii^ I, 17. 

(A2) if Xi = 0, then fr = 0; in particular, if x e X,, then = for / = 1 , 9. 
(A3) = if / > 9; there are no new infections in uninfected compartments. 
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(A4) if X e X,, then = and f^ix) = for / = 1 , 9. 

(A5) if is set to 0, then all eigenvalues of Df(xo) have negative real parts. 

To construct the next generation matrix, we only consider infected and exposed compartments. The equations are trans- 
formed as follows. 



d 
dt 



/321S1I2/N2 




I323S3I2/N2 


PuSaIxINx +.fl32ASAl2lN2 +/334S4I3/N3 




diEiNi/Ki +E1E1 
-exQi+dxhNilKx-siEi 

d2E2N2/K2 + S2E2 
-£2 £2 + d2hN2lK2 + 72^2 + y"2/2 

diEjNi/Ki + e-^Et, 
d-ihNilK-i- SiE-i 
dAE^NiilKA + S4E4 

d4l4N4/K4 — S4E4 + 74/4 + JU4/4 



The equation system is nonlinear; we linearize it, deriving the two Jacobian matrices. First, the partial derivative of ^ with 
respect to each variable at the disease free equilibrium is as follows [16]. 










/3i2^ 



biqi 















A4^ 

' 





y82i^ 
' 



s" 

/323^ 

' 

s" 

//?24^ 










y632 





p 







rO 

34 73) 



























(74) 



Second, the partial derivative of y with respect to each variable at disease free equilibrium is as follows. 



V : 
































bi +si 























-Oi 


-si 


bi 





























b2 + B2 


























-S2 


^2 + 72 + ^J■2 





























bi +£3 


























-£3 


bi 





























b4 + E4 


























-E4 


^4 + 74 + yU4 



The inverse of matrix V is computed as foUows. 





X 












1 

ei 














_1_ 















1 

bi+e2 

£2 

(b2+S2)(hi+y2+H2) 








02 +72+^12 













1 

i>3+£l 

£3 

by{bi+e^) 














X 

A3 
















































1 

/74+£4 






g4 



1 



(i4+£4)(fe4+y4+/i4) bi+Ji+lH 



(75) 



Finally, the next generation matrix, which is the product FV is as follows. 
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r A 1 


R, 


A , 
^1 

















' 


U 


U 


U 




L>1 


U 


U 


U 


u 





























El 


Fx 


Gi 








Hi 


h 












































Jl 


Li 









































Ml 




^1 


Gi 


^1 


5i 






































where: 



Ai = qi 

Ci 
Di 
El 
Fi 
Gi 
Hi 
h 

Jl 
Li 
Ml 

Ni 

Pi 

Qi 
Ri 
Si 

Ti 



bi +si 

S2fi2lblKid2 

(bi + siXbi +72+ l^2)dib2K2 
fi2ibiKid2 

0)2+72 +H2)dlb2K2 
Pl2b2K2di 

b\d2Ki 
si/3i2b2K2di 
(bi +Ei)b\d2Ki 
/3i2b2K2di 
b\d2Ki 

g3/?32fc2^2'i3 

{bi + sj,)bld2Ki 

^32b2K2dj, 

bld2K, 

(b2 + S2)(^2 +72+ l^2)dib2K2 

bsKjd2fi23 
0)2 +72+ H2)dib2K2 

biiK^diPn 
b\dAK\ 
SiPinb^Kudi 



(bi + si)b]d4Ki 

b^diKiiPn 
b\dAKi 

fS2b4K4d2/324 

ib2 + S2){b2 +72+ H2)d4b2K2 

fb4K4d2fi24 
0)2 + 72+ f^2)d4b2K2 
S3b4K4d3/3s4 



(bs + Si)bld4Ki 

b4K4diPi4 
blK3d4 



Recall that the reproduction number is the spectral radius of FV ^ we compute the nine eigenvalues A, of FV ^ to select the one 
with maximum magnitude. 
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IFV"' - AI\ = 



A\ — A 


oi 


A 1 























-A 





r, 
^ 1 






















-A 




















£i 


^'i 


Gi 







H, 


/i 




















-A 

























Li 


-/I 





































Ml 




Pi 


fil 


Ri 


5i 


Ti 


-A 





























-A 



= -A^[A^ -AiA^-(CiFi + JiHi)A + (HiAiJi +A1C1F1 -BiCi^i)] (76) 



B,C,E, = qiflfM^ ^77^ 

bi(bi + ei)(b2 + £2X^2 +1^2 + 72) 

A,C,F, = il£if«^ (78) 

+ ei)(/72 + £2X^2 +A'2 + 72) 

Because BiCiEi = AjCiFi as shown in Equation (78) and (77) , equation (76) can be rewritten as follows. 



-A^[A^ - AiA^ - (CiFi + JiHi)A + H^AiJi] = (79) 
Equation (79) has six zero roots. We only need to solve the following equation to find max (/ =1,2, 3). 



A^ -AiA^ -(CiFi +JiHi)A + HiAiJi =0 (80) 



Equivalently, 



^3 _ j^2_j SiS2/32l/3l2 _^ £253^23^32 

blibi + Si)ib2 + £2X^2 +f^2+ 72) ^3(^3 + £3X^2 + £2)(^2 +1^2 + 72) 

?l£2£3^23^32 (gj) 



b^ib^ + Si)(b2 + £2)(^2 +1^2+ 72) 



We calculate the reproduction number numerically with 5000 different sets of parameters uniformly distributed within the range 
in [18]. The histogram of the reproduction number is shown in Figure 7. From the histogram, we can see that Rq can be greater 
or smaller than 1. In particular, the mean is 1.17 and the maximum is 3.68, respectively. 




12 3 4 

Possible values of RO 



Figure 7: Histogram of the reproduction number, tlie mean is 1.17, the maximum is 3.68 
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Upper and Lower Bound for Rq 

Although we are able to only obtain the exact expression of Ro numerically, we determine the lower bound and the upper 
bound of Ro in the following. 



\FV-^ - AI\ 
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Where: 
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Jl 
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<?1E1 








Ai-^ Bi 
-/I Ci 







El 






S2l32lhiK]d2 

+E2 )(''2 +72 +/'2 )dl ''2 ^2 



£2^3^3^2)323 
(^2 +E2 )(^2 +72 +/^2)'^3 ^2 ^^2 



-.1 i/l 






(b^+ej,)b\d2K2 





-a^\:r-ai\ 



Matrix FV ' has five zero roots. To find out the spectral radius of FV we only need to find out the spectral radius of matrix 
J[ which can be rewritten as follows. 






P\2b2K2di 

h]d2Ki 








Sil3\2h2K2d\ 
{bi+eOb]d2Ki 







£2P2\b\K\d2 



(b2+S2)(b2+y2+P2)dtb2K2 








£3/^32^2 A^2'j3 

{hi+Bi)b\d2Ki 





q\ 







*1+£| 














(i2+£2)(^'2+y2+/J2K3''2^^2 

S is the first matrix and C is the second matrix. The matrix C has three multiple eigenvalues A2 - At, - Aa, — 0. Since 
Rank (Ail - C)|,i,=^,=^j=o - 1, there are three linear independent eigenvectors corresponding to zero eigenvalue. Overall, matrix 
C has four linear independent eigenvectors. Therefore, matrix C can be diagonalized as follows. 
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where: 
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Substitute C in the expression of = S + C. 
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)P-^ = P{X + Y)P-^ 



where 






P2\b2K2d\ 
b\d2Ki 





B,e2fi2lb,Kid, 



(A 1 +ei )(fe2 +£2 Xfc +72 +H2)d\ b2 K2 
e2/32ibiKid2 



q\ib2+E2)(b2+y2+^2)d\b2K2 



S2P2ibiK^d2 






Eil3i2h2K2di 



(bi+E2)d2b3^K, 





(b2+B2){b2+-y2+f2)d?.b2K2 

Matrix 7 is a 4 x 4 nonnegative diagonal matrix [Yjj = for / j) and ^ F,, < qi 










< 00, (/ = 2, 3,4), and matrix 



X is nonnegative. Therefore, p(X) ^ p(X + Y) ^ p{X) + qi according to Theorem 1. in [10]. The eigenvalues of matrix X + Y 
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are the same as those of J{ because the two matrices are similar. Similarly, the eigenvalues of matrix X are the same as those of 
matrix S . 



p{X + Y) 
P(X) 



Ro 



S2 
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\ (^2 +£2X^2 +72 +jU2)'^^i(^i ^73(^3 +£3) 



If we only count the horizontal transmission and denote the new F (resp. V) by Fh (resp. Vh) , Fh and Vh are as follows. 



By calculation, 



V„ = 
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^1 (b2 + E2)(b2 + 72 + mVbiibi +si) biibs+Si) 



Therefore, 



< s$ /?o < < + qi (82) 
We denote and R^ + qi as R^ and R^, respectively. 



Numerical Comparison among Rq, and R^ 

Rq and Rq + q\ are the lower bound and the upper bound of Rq, respectively. Therefore, R^ + q\ < \ ^ Rq < \, and 
R^^ > 1 ^ > 1 . To verify that the derived bounds are tight, we perform extensive simulations using 5000 sets of parameters 
uniformly distributed within the range defined in [18]. In Figure 8(a) and 8(b), R\!^ - Rq + qi vs. Ro and 7?^ = R^ vs. Rq are 
plotted, respectively. First, the difference between the exact values and each bound is very small. In fact, the red and green points 
lay very close to the line y - x. Additionally, in the case of the upper bound, the red points are just slightly above the line y - x, 
while in the case of the lower bound, the green points are just slightly below the line y - x. 
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(b) The reproduction number and its lower bound. 
Figure 8: The reproduction number and its upper and lower bound 
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Biological Interpretation of Bounds for Rq 

The bounds for Rq, as given in inequalities (26), can be interpreted biologically as follows. The lower bound, R^, is the 
reproduction number for horizontal transmission because R^ - p (F/zV^'), where p (F/zV^') represents the spectral radius of 
the next generation matrix for horizontal transmission FhVJj^ . The upper bound is given by the sum of R^ and a second term 
that is only related to vertical transmission, i.e. from mothers to their offspring in the Aedes mosquito population. 

R^ includes Aet/ei-livestock interaction and CMfex-livestock interaction. More specifically, we define the reproduction num- 
ber due to the interaction between Aedes and livestock represented by R^^^^^^ as 



H(A-L) 




Si 



\(b2 + S2)(b2+y2+l^2)bi(bi +Si) 



Rq^^ ^' can be rewritten as follows. 



- i^^^l^F^ ^] (83) 

\ dx-^ {dx-^ + sx) (d2^ +y2+fJ-2)(d2Y; +f^2) 



because 



bx^dx — 
Kx 

bi - d2 — 

K2 



where: 

A^* - the total number of Aedes mosquitoes at disease free equilibrium. 
A^2 - the total number of livestock at disease free equilibrium. 

Rn^^^^^ consists of the product of four terms. Each infected Aet/ei mosquito can infect susceptible livestock throughout 
its lifetime. Similarly, each infected livestock can infect — susceptible Aedes mosquitoes during its lifetime. The prob- 

ability of Aedes mosquitoes and livestock surviving through the incubation period to the point where they become infectious 
is — and — ^ — , respectively. Therefore, i?^^'*"^' is the geometric mean of the average number of secondary livestock 

infections produced by one Aedes mosquito vector in the first square bracket in (83), and the average number of secondary Aedes 
mosquito vector infections produced by one livestock host in the second square bracket in (83). 

Similarly, we define the reproduction number due to the interaction between Culex and livestock represented by iJ?^''"^^ as 



pH(C-L) 
^0 



B2 63^632^: 
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^ (^2 + e2)(^2 +72+ f^2) biibi + £3) 



We can rewrite R^^"^ ^' as follows. 



dH(c-l) _ Ir 1^32 fi3 ir g2 1 

\ di ("3 ^ + £3) (122 + 72 + /"2) (02 + £2) 

because 

m 

b,=d3-^ 

where: 

- the total number of Culex mosquitoes at disease free equihbrium. 
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R 



H{C-L) 




also consists of the product of four terms. Each infected Culex mosquito can infect susceptible livestock 



throughout its lifetime. Similarly, each infected livestock can infect ^. " susceptible Culex mosquitoes. The probability of 

Culex mosquitoes surviving through the incubation period to the point where they become infectious is — — . Similarly, the 

Therefore, 



probability of livestock surviving through the incubation period to the point where they become infectious is 



R 



H(C-L) 




is the geometric mean of the average number of the secondary livestock infections produced by one Culex mosquito 
vector in the first square bracket in (84), and the average number of secondary Culex mosquito vector infections produced by 
one livestock in the second square bracket in (84). 

The expression (27) for , can be rewritten as R^ = ■yJiR"''^'^^)^ + (^^ *'^"'^')^ ' where the dependence of R^ on 
and R^^^^^^ is shown in Figure 9. The square root is due to the vector-host-vector viral transmission path [18? , 33]. Obviously, 
the horizontal reproduction number increases with the increase of each of the four terms in R^^^^^^ and r'^^^^^\ 




H(A-L) 



Figure 9: The interpretation of 
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